#define PERFORMANCE_TESTING  0

//
//  Standard headers
// 
#include <iterator>   // for stream_iterators, inserter
#include <iostream>
	using std::cout;
	using std::cerr;
	using std::cin;
	using std::ofstream;
	using std::ifstream;
#include <string>
	using std::string;
#include <vector>
	using std::vector;
#include <map>
	using std::map;
#include <utility>
#include <algorithm>
	using std::min;
#include <numeric>
#include <stdexcept>
	using std::runtime_error;
#include <boost/range/algorithm_ext/push_back.hpp>
    using boost::push_back;
#include <limits>
	using std::numeric_limits;
#include <algorithm>
	using std::max;

//
//  My common headers
// 
#include <progress_indicator.h>
#include <commify.h>

//
//  Boost headers
// 
#include <boost/algorithm/string/join.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>
#include <boost/tr1/tuple.hpp>
#include <boost/tr1/unordered_map.hpp>
	using std::tr1::unordered_map;
#include <boost/unordered_set.hpp>
	using boost::unordered_set;
#include <boost/range/algorithm.hpp>
	using boost::sort;
#include <boost/format.hpp>
		using boost::format;
		using boost::io::group;
#include <boost/algorithm/string/replace.hpp>
		using boost::algorithm::replace_all;
#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
		namespace io = boost::iostreams;

#include <boost/assign/list_of.hpp> // for 'list_of()'
		using namespace boost::assign; // bring 'list_of()' into scope

#include <boost/assert.hpp> 
#include <boost/lexical_cast.hpp> 
#include <boost/spirit/include/karma.hpp>
#include <boost/regex.hpp>

#include <boost/random/mersenne_twister.hpp>
	using boost::mt19937;
#include <boost/random/uniform_int.hpp>
	using boost::uniform_int;
#include <boost/random/variate_generator.hpp>	
	using boost::variate_generator;




//
//  tr1 headers
// 
#include <boost/tr1/tuple.hpp>
	using std::tr1::tuple;
	using std::tr1::make_tuple;
	using std::tr1::get;
	
#include "tbb/task_scheduler_init.h"
    // tbb::task_scheduler_init

#include "intervals_per_contig.h"
	using intervals_per_contig::t_interval;
	using intervals_per_contig::t_intervals;
	using intervals_per_contig::t_intervals_per_contig;
	using intervals_per_contig::t_vec_intervals;
	using intervals_per_contig::t_map_str_intervals;
#include "read_options.h"

#include "tbb_processor.h"











//!_______________________________________________________________________________________
//!     
//!		get_intervals_from_all_files
//!			converts from intervals per contig to genomic coordinates
//!     
//!		\author lg (1/11/2011)
//!     
//!     \param opt 
//!     \param intervals 
//!     \param contig_genome_offsets 
//!_______________________________________________________________________________________
void get_intervals_from_all_files(prog_options& opt, vector<t_intervals>& intervals, t_contig_genome_offsets& contig_genome_offsets)
{
	for (unsigned i = 0; i < opt.interval_file_names.size(); ++i)
	{

		//
		//	Read intervals for each file
		// 
		string file_name = opt.interval_file_names[i];
		opt.both() << "Reading intervals from " << file_name << "...\n";
		io::stream<io::mapped_file_source> f_intervals(file_name);
		if (!f_intervals)
			throw std::logic_error("Can't open " + file_name);
		t_intervals_per_contig temp_intervals_per_contig;
		read_intervals(	f_intervals, temp_intervals_per_contig, opt.zero_based_intervals /*one_based_coordinates?*/);

		// 
		// 	Convert	from per contig coordinates to per genome coordinates
		intervals.push_back(t_intervals());
		t_vec_intervals& curr = intervals.back().data;
		BOOST_FOREACH(t_map_str_intervals::value_type& k_v, temp_intervals_per_contig.data) 
		{
			long long contig_pos_to_genome = contig_genome_offsets[k_v.first];
			BOOST_FOREACH(t_interval& interval, k_v.second.data) 
				curr.push_back(t_interval(	interval.beg + contig_pos_to_genome,
											interval.end + contig_pos_to_genome));
		}


		opt.both() << "    " << commify(curr.size()) << " intervals\n";
		intervals.back().combine_overlapping();
		opt.both() << "    " << commify(curr.size()) << " non-overlapping intervals\n";
		opt.both() << "    " << commify(intervals.back().len()) << " bp\n";
	}
	opt.verbose(1) << "Finished reading intervals...\n";

}

//!_______________________________________________________________________________________
//!     
//!		get_offsets_into_genome_pos_per_contig
//! 		Saves the offsets required per contig to convert contig coordinates to
//! 		genomic ones.
//!     
//!		\author lg (1/11/2011)
//!     
//!     \param opt 
//!     \param contig_genome_offsets 
//!     
//!     \return int total genome size
//!_______________________________________________________________________________________

struct is_canonical_contig
{
	bool operator()(const string& contig)
	{
		static boost::regex canonical_only("([0-9]+)|X|x");
		return regex_match(contig, canonical_only);
	}
};
long long get_offsets_into_genome_pos_per_contig(prog_options& opt, t_contig_genome_offsets& contig_genome_offsets)
{
	unordered_map<string, t_interval> contig_sizes;

	string file_name = opt.contig_sizes_file_name;
	opt.both() << "Reading contig sizes  from " << file_name << "...\n";

	io::stream<io::mapped_file_source> f_intervals(file_name);
	if (!f_intervals)
		throw std::logic_error("Can't open " + file_name);
	t_intervals_per_contig sizes_per_contig;
	read_intervals(	f_intervals, sizes_per_contig, opt.zero_based_intervals /*one_based_coordinates?*/);
	vector<string> contig_names = sizes_per_contig.get_sorted_contigs();


	unsigned cnt_sizes = sizes_per_contig.size();
	if (cnt_sizes != contig_names.size())
		opt.both() << "WARNING! Some contigs have more than one size\n";

	//	only keep 1-99,X
	//contig_names.erase(std::partition(
	//								contig_names.begin(), contig_names.end(), is_canonical_contig()), 
	//					contig_names.end());
	opt.both() << "    " << contig_names.size() << " contigs\n";
	opt.log() << "    contigs = [ " << boost::join(contig_names, ", ") << " ]\n";


	// only take the first size
	long long genomic_offset = 0;
	BOOST_FOREACH(const string& contig_name, contig_names) 
	{
		t_interval interval = sizes_per_contig.data[contig_name].data[0];
		genomic_offset -= interval.beg;
		contig_genome_offsets[contig_name] = genomic_offset;
		//opt.verbose() << boost::format("contig =%1$14s, beg = %2$15s, end  = %3$15s, offset = %4$15s\n") 
		//								% contig_name
		//								% commify(interval.beg)
		//								% commify(interval.end)
		//								% commify(genomic_offset);
		genomic_offset += interval.end;
	}
	return genomic_offset;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// 
//      Main
// 
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888 
//!_______________________________________________________________________________________
//!     
//!		output_counts
//!     
//!		\author lg (3/3/2011)
//!     
//!     \param opt 
//!     \param counts_per_category_per_interval_type 
//!     \param category_names 
//!_______________________________________________________________________________________
void output_counts(prog_options& opt, const std::vector<std::vector<unsigned> >& counts_per_category_per_interval_type, vector<string>& category_names)
{
	if (opt.get_verbosity() == 0)
		return;

	// 
	// get format strings for categories
	// 
	vector<string> category_format_strs;
	// size of string with commas every thousand
	string::size_type min_number_len = 12;
	for (unsigned ii = 0; ii < category_names.size(); ++ii)
	{
		unsigned category_name_len = std::max(min_number_len, category_names[ii].length()) + 1;
		category_format_strs.push_back((format("%%1$%1$d") % category_name_len).str());
	}

	// 
	// Get format strings for categories so each is the correct width
	// 
	std::size_t interval_name_len = min_number_len;
	for (unsigned ii = 0; ii < opt.interval_type_names.size(); ++ii)
		interval_name_len = max(interval_name_len, opt.interval_type_names[ii].length());
	string interval_format_str = (format("%%1$%1$ds") % (interval_name_len + 1)).str();



	//
	//  Print counts
	// 

	//  Category names
	cout		<< format(interval_format_str) %  "";
	for (unsigned ii = 0; ii < category_names.size(); ++ii)
		cout		<< format(' ' + category_format_strs[ii] + 's') % category_names[ii];
	cout		<< std::endl;

	//  Numbers and accumulate total
	vector<unsigned> totals(category_names.size(), 0U);
	for (unsigned ii = 0; ii < opt.interval_type_names.size(); ++ii)
	{
		cout		<< format(interval_format_str) %  opt.interval_type_names[ii];
		for (unsigned jj = 0; jj < category_names.size(); ++jj)
		{
			cout 		<< format(' ' + category_format_strs[jj] + 's') % commify(counts_per_category_per_interval_type[jj][ii]);
			totals[jj] += counts_per_category_per_interval_type[jj][ii];
		}
		cout << "\n";
	}

	// 
	// Print Totals
	// 
	cout		<< format(interval_format_str) %  "";
	for (unsigned ii = 0; ii < totals.size(); ++ii)
		cout 		<< format(' ' + category_format_strs[ii] + 's') % string(min_number_len, '-');
	cout << "\n";
	cout		<< format(interval_format_str) %  "Total";
	for (unsigned ii = 0; ii < totals.size(); ++ii)
		cout 		<< format(' ' + category_format_strs[ii] + 's') % commify(totals[ii]);
	cout << "\n\n";
	cout << std::endl;


	// 
	// Proportions
	// 
	//  Category names
	cout << "Proportions:\n";
	cout		<< format(interval_format_str) %  "";
	for (unsigned ii = 0; ii < category_names.size(); ++ii)
		cout		<< format(' ' + category_format_strs[ii] + 's') % category_names[ii];
	cout		<< std::endl;

	//  Numbers and accumulate total
	for (unsigned ii = 0; ii < opt.interval_type_names.size(); ++ii)
	{
		cout		<< format(interval_format_str) %  opt.interval_type_names[ii];
		for (unsigned jj = 0; jj < category_names.size(); ++jj)
			cout 		<< format(category_format_strs[jj] + ".2f%%") % 
									(counts_per_category_per_interval_type[jj][ii] * 100.0 / totals[jj]);
		cout << "\n";
	}
	cout << "\n";

	cout << std::endl;

}

//!_______________________________________________________________________________________
//!     
//!		output_counts_to_tab_delimited_summary
//!     
//!		\author lg (3/3/2011)
//!     
//!     \param opt 
//!     \param counts_per_category_per_interval_type 
//!     \param category_names 
//!_______________________________________________________________________________________
void output_counts_to_tab_delimited_summary(prog_options& opt, const std::vector<std::vector<unsigned> >& counts_per_category_per_interval_type, vector<string>& category_names)
{
	for (unsigned ii = 0; ii < category_names.size(); ++ii)
		opt.summary() << "\t" << category_names[ii];
	opt.summary() << "\n";
	vector<unsigned> totals(category_names.size(), 0U);
	for (unsigned ii = 0; ii < opt.interval_type_names.size(); ++ii)
	{
		opt.summary() << opt.interval_type_names[ii];
		for (unsigned jj = 0; jj < category_names.size(); ++jj)
		{
			opt.summary() << "\t" << counts_per_category_per_interval_type[jj][ii];
			totals[jj] += counts_per_category_per_interval_type[jj][ii];
		}
		opt.summary() << "\n";
	}
	// 
	// Totals
	// 
	opt.summary() << "Total";
	for (unsigned ii = 0; ii < totals.size(); ++ii)
		opt.summary() << "\t" << totals[ii];
	opt.summary() << "\n";
	opt.summary() << "\n";
	opt.summary() << std::endl;

	// 
	// Proportions
	// 
	for (unsigned ii = 0; ii < category_names.size(); ++ii)
		opt.summary() << "\t" << category_names[ii];
	opt.summary() << "\n";
	for (unsigned ii = 0; ii < opt.interval_type_names.size(); ++ii)
	{
		opt.summary() << opt.interval_type_names[ii];
		for (unsigned jj = 0; jj < category_names.size(); ++jj)
		{
			opt.summary() << "\t" << counts_per_category_per_interval_type[jj][ii] * 100.0 / totals[jj] << "%";
		}
		opt.summary() << "\n";
	}
	opt.summary() << "\n";
	opt.summary() << std::endl;
}


int main (int argc, char *argv[])
{
	try
	{
		// 
		//  Read options
		//  	return if user chose help or version
		// 
		prog_options opt;
		if (!read_options(opt, argc, argv))
			return 0;




		//
		//	Get contig sizes 
		// 
		t_contig_genome_offsets contig_genome_offsets;
		unsigned genome_size = get_offsets_into_genome_pos_per_contig(opt, contig_genome_offsets);



		//
		//  Read intervals from all files
		// 
		vector<t_intervals> intervals;
		get_intervals_from_all_files(opt, intervals, contig_genome_offsets);

		// 
		// add whole genome to intervals as catchall
		//
		intervals.push_back(t_intervals());
		intervals.back().data.push_back(t_interval(0, numeric_limits<unsigned>::max()));
		opt.interval_type_names.push_back(opt.default_interval_type_name);


		//
		//  start tbb
		// 
        tbb::task_scheduler_init init( opt.jobs );


		// 
		//	both files 
		// 
		
		ifstream input_stream;
		//io::stream<io::mapped_file_source> input_stream;
		
		input_stream.open(opt.input_file_name.c_str());
		if (!input_stream)
			throw std::logic_error("Can't open " + opt.input_file_name);
		opt.both() << "Input= " << opt.input_file_name << std::endl;

			
		// 
		// Parse all per genotype data per strain
		// 
		opt.verbose() << "Parsing lines..." << std::endl;
		unsigned cnt_threads_parallel_tasks = opt.jobs * 2;
		vector<vector<unsigned> >				counts_per_category_per_interval_type;
		vector<string>							category_names;
		unsigned cnt_lines = parse_lines_pipeline(	cnt_threads_parallel_tasks, 
													input_stream, 
													opt.verbose(),
													contig_genome_offsets,	
													opt.zero_based_snps,
													opt.intervals_only,
													counts_per_category_per_interval_type,
													category_names,
													intervals);

		//
		//	sort counts by category name
		// 
		map<string, vector<unsigned> > map_names_to_counts;
		for (unsigned ii = 0; ii < category_names.size(); ++ii)
			map_names_to_counts[category_names[ii]] = counts_per_category_per_interval_type[ii];
		sort(category_names);
		for (unsigned ii = 0; ii < category_names.size(); ++ii)
			counts_per_category_per_interval_type[ii] = map_names_to_counts[category_names[ii]];

		//
		//	output counts 
		// 
		output_counts_to_tab_delimited_summary(opt, counts_per_category_per_interval_type, category_names);
		output_counts(opt, counts_per_category_per_interval_type, category_names);

		//
		//  Finished parsing
		// 
		opt.both() << commify(cnt_lines) << " lines parsed.\n\n";


	}
	catch(const std::exception& e)
	{
		// indent error string
		string error_string = e.what();
		replace_all(error_string, "\n", "\n\t");
		std::cerr
				<< "\n"<< string(60, '=') << "\n"
				<< "\nError:\n\t" << error_string 
				<< "\n\n" << string(60, '=') << "\n\n";
		return -1;
	}
	return 0;
}


